#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <unistd.h>
#include <assert.h>
#include <sys/time.h>
#include <omp.h>

//#define N 258
double A[N][N];
#ifdef CHECK
double test[N][N];
#endif
int bs;

int max(int a, int b)
{
	if(a>b) return a;
	return b;
}

int min(int a,int b)
{
	if(a<b) return a;
	return b;
}   

int gauss_seidel(const int bs,const int L)
{
	int i,j,ii,jj;
	int iter;
   
    for(iter = 0; iter < L; ++iter)
    {
        for(ii = 1; ii <= N -2; ii += bs)
        {
            for(jj = 1; jj <= N - 2; jj += bs) 
            {
                //need omp4.0
                #pragma omp task depend(inout:A[ii:min(N-1-ii,bs)][jj:min(N-1-jj,bs)]) depend(in:\
				A[max(ii-bs,0):(ii-bs>0?bs:1)][jj:min(bs,N-1-jj)],\
				A[ii:min(bs,N-1-ii)][max(jj-bs,0):(jj-bs>0?bs:1)],\
				A[ii+min(bs,N-1-ii):min((N-1-ii-bs>0?N-1-ii-bs:1),bs)][jj:min(bs,N-1-jj)],\
				A[ii:min(bs,N-1-ii)][jj+min(bs,N-1-jj):min((N-1-jj-bs>0?N-1-jj-bs:1),bs)])
				for(i=ii;i<ii+(N-1-ii>bs?bs:N-1-ii);i++) 
						for(j=jj;j<jj+(N-1-jj>bs?bs:N-1-jj);j++)
						{
							A[i][j] = 0.2*(A[i][j] + A[i-1][j] + A[i+1][j] + A[i][j-1] + A[i][j+1]);
						}
            }
        }
        
    }
    return 0;
}



void Init()
{
		int i,j;
		for(i = 0;i < N;i++)
			for(j = 0;j < N;j++)
			{
				A[i][j] = (double)(i*2.5 + j*3.3 ) /3.0;
			
			}
				
}
#ifdef CHECK
int check(const int ITER,const int bs)
{
	int i,j,itr; 
	printf("begin check\n");
    for(i = 0; i < N; i++) 
    for(j = 0; j < N; j++)
          test[i][j] = (double)(i*2.5 + j*3.3 ) /3.0; 
   for(itr = 0; itr < ITER; ++itr)
    {
        for(i = 1; i < N-1 ; i++)
        for(j = 1; j < N-1 ; j++)
           	test[i][j] = 0.2* (test[i][j] + test[i-1][j]+test[i+1][j] + test[i][j-1] + test[i][j+1]);
    }

    for(i = 1; i < N - 1; i++)
    for(j = 1; j < N - 1; j++)
       if( fabs(test[i][j]  - A[i][j]) > 1e-11 )
       {
		   printf("i : %d,j : %d\n",i,j);
		   printf("test : %lf\n",test[i][j]);
		   printf("A : %lf\n",A[i][j]);
           return 1;
       }
	printf("correct\n");
    return 0;
}
#endif

int main(int argc, char** argv)
{
    int num_threads;
    struct timeval start;
    struct timeval end;
    double total_time;       
    
    if(argc > 4)
    {
        printf("input arg error.Usage: ./a bs iter num_threads\n");
        return 0;
    }
    num_threads = 1;
    bs = 16;
    int iter = 1;
    if(argc > 1)
        bs = atoi(argv[1]);
    if(argc > 2)
        iter = atoi(argv[2]);
    if(argc > 3)
        num_threads = atoi(argv[3]);
    omp_set_num_threads(num_threads);
    Init();
    int max_threads = omp_get_max_threads(); 
	gettimeofday(&start,NULL);
    #pragma omp parallel num_threads(max_threads)
	#pragma omp single
    gauss_seidel(bs,iter);
    gettimeofday(&end,NULL);   
     
    total_time = 0.0;
    total_time+=(end.tv_sec-start.tv_sec)+(end.tv_usec-start.tv_usec)/1000000.0;
#ifdef CHECK
    
    if(check(iter,bs))
    { 
        printf("check_error\n");
        exit(1);
    }  
#endif
    printf("total time: %lf\n",total_time);
    return 0;
}
